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Abstract 

In this paper we show that a methodology based on a sampling 
with the Gaussian function of kind / {c^/n), where c and h 

are some constants, leads to the Fourier transform that can be rep¬ 
resented as a weighted sum of the complex error functions. Due to 
remarkable property of the complex error function, the Fourier trans¬ 
form based on the weighted sum can be significantly simplified and 
expressed in terms of a damping harmonic series. In contrast to the 
conventional discrete Fourier transform, this methodology results in a 
non-periodic wavelet approximation. Consequently, the proposed ap¬ 
proach may be useful and convenient in algorithmic implementation. 

Keywords: complex error function, Faddeeva function, Fourier trans¬ 
form, sampling, Gaussian function, numerical integration 


1 Introduction 

Sampling is an efficient methodology that can be used for functional approx¬ 
imations DEI- It is based on interpolation through a set of discrete points 
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{tu) f (^n)} selected from a function / (t) within an interval [a, b]. The sine 
function 

{ sinf 
t ’ 

is one of the most popular in sampling [31 HE]. Specihcally, due to its 
remarkable sampling property / {t) can be represented as a series approxi¬ 
mation HE] 

N 

f{t)- X] sine /(tn), (1) 

n=—N 

where h is a small parameter that at = nh can be regarded as the step. 
We have shown recently a new methodology of sampling based on incomplete 
cosine expansion of the sine function in numerical integration [6]. 

Another very useful sampling tool can be constructed by using the Gaus¬ 
sian function . The sampling property of the Gaussian function can 
be observed by changing its variable as t —)■ t/e and then considering the 
limit 

hme“^*/^^ = 5 (t), 

£—^0 

where the Kronecker’s delta 5 {t) is a building block of the shah function, also 
known as the comb function: 


t ^ 0 
t = 0, 


OO 

Ul{t)= Y, Ht-tn) 

n=—oo 


(see application of the shah function III (f) in sampling in [T]). Gonsequently, 
this example shows that the Gaussian function can also be utilized for func¬ 
tional approximations by sampling methodology. 

It is convenient to modify the Gaussian function by multiplying a 
constant h/ {cy/Tr) and making change of the variable as t ^ t/c. This results 
in 



( 2 ) 


where the c is the htting parameter. Figure 1 shows three variations of the 
sampling Gaussian function (|^ computed at h = 0.25, c = 0.15 (black curve), 
h = 0.25, c = 0.2 (blue curve) and h = 0.25, c = 0.25 (red curve). We may 
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attempt to approximate / (t) by replacing the sine functions in equation Q 
with the Gaussian functions of kind ([^ as given by equation 


fit) 


h 


N 


c^/tt ^ 

* n — 




fitn 


n=-N 


Consequently, choosing the sampling points separated equidistantly with step 
h between two adjacent points we obtain 

h ^ 

/(i)«-G ^ (3) 

C-v/TT 

^ n=-N 


The Gaussian functions 



Fig. 1. The Gaussian sampling functions computed at h = 0.25, 
c = 0.15 (black curve), h = 0.25, c = 0.2 (blue curve) and h = 0.25, 
c = 0.25 (red curve). 


The application of the Gaussian function of kind (|^ can be justihed by 
considering a simplest case when / (t) = 1. Figure 2 shows the corresponding 
approximations at h = 0.25, c = 0.15, iV = 10 (blue curve) and h = 0.25, c = 
0.2, = 10 (red curve). The multiplication of the number of the sampling 

points to the step determines the range for approximation {2N + 1) h = 5.25. 
As a consequence, the approximated functions shown in Fig. 2 resemble the 
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window function. In particular, these functions are approximately equal to 
1 within the range t G [—5.25/2,5.25/2] and to 0, otherwise. As we can 
see from this hgure, at smaller c the approximated function oscillates near 
the constant 1 at the top of the curve stronger (blue curve). However, as c 
increases the oscillation rapidly decreases (red curve) and practically vanishes 
at c ^ h. This signihes that if the change in curve between any two adjacent 
sampling points is relatively small, then any function can be approximated 
by sampling with the Gaussian function ([^. 


Window function approximations 



Fig. 2. The window function approximations computed at h = 0.25, 
c = 0.15, = 10 (blue curve) and h = 0.25, c = 0.2, = 10 (red 

curve). 


In this work we introduce an application of the complex error function to 
the Fourier analysis. In particular, we show that the use of equation (|^ in the 
Fourier integration leads to a weighted sum of the complex error functions. 
Due to remarkable property of the complex error function this approach 
provides efficient computational methodology in the Fourier transform as a 
damping harmonic series. 
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2 The complex error function 


The complex error function, also known as the Faddeeva function or the 
Kramp function, is defined as [3 El El ED] 


2: 



0 


where z = x + iy is the complex argument. This function is a solution of the 
following differential equation [10] 


w' {z) + 2zw {z) 



m;(0) = 1. 


The complex error function hnds broad applications in many helds of Ap¬ 
plied Mathematics [H [9l E], Physics [1211131 Ell ESI ESI EH EB] and Astron¬ 
omy [19]. In Applied Mathematics it is closely related to the error function 
of complex argument ra 

w (iz) = e~^'^ [1 — erf {—iz)] erf {z) = 1 — e~^^w (iz ), (5) 


the normal distribution function m 


$(z) = 






1 _ 


the Fresnel integral PD] 


F(2) = / 

0 

= (1 + z) [1 - (1 + i) z/2)\ /2 

and the Dawson’s integral [H El El El El 121] 

—e~^ + w {z) 

2i ■ 


daw (z) = e ^ du = \/n 
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In Physics and Astronomy the complex error fnnction is related to the 
Voigt fnnction [ini [121 [131 [El [El IE] 





that describes the spectral behavior of the photon emitting or absorbing 
objects (photo-lnminescent materials [E] , planetary atmosphere [El El El 
El E], celestial bodies [E] and so on). Specihcally, the Voigt function 
represents the real part of the complex error function 

K (x, y) = Re [w (x, y)], y> 0. 

Other functions that can be expressed in terms of the complex error function 
are the plasma dispersion function m, the Gordeyev’s integral [E], the 
rocket flight function ESI and the probability integral |2nj . 

The complex error function Q can be represented alternatively as (see 
equation (3) in [23] and [2^) see also Appendix A in [22] for derivation) 


OO 



0 


Using the change of the variable as m —)■ 2m in the integral above leads to [28] 



0 


Further, we will use this equation in derivation of the weighted sum. 

A rapid C/C++ implementation (RooFit package from CENR’s library) 
for computation of the complex error function with average accuracy ~ 10“^® 
has been reported in the recent work [29] . 
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3 Derivations 


3.1 Weighted sum 

There are several definitions for the Fourier transform. In this work we will 
use the following definitions [1] 

OO 

.^{/(«)}M= j (7) 

— OO 


and 


^ {F{y)]{t)= / F (i/) 


( 8 ) 


Thus, the relationships between functions / (t) and F (z/) are performed by 
two operators F {/ (t)} (z/) = F (z/) and F~^ {F (i/)} (t) = f (t) correspond¬ 
ing to the forward and inverse Fourier transforms, respectively. In signal 
processing the arguments t and z/ in these reciprocally Fourier transformable 
functions / {t) and F (z/) are interpreted, accordingly, as time vs. frequency 
0. 

We can find an approximation to the Fourier transform of the function 
/ (t) by substituting approximation (|^ into equation ([^. This leads to 


FM 


h 


C0f 


h 


' N 

E' 

n=-N 


' t — nh \ ^ 


) f{nh). 


—27riiyt 


dt 


Cx/TT 


N ^ 

E 

_ 


^ 9 nht / n h \ 2 

e ^ e 


) /(nh) 




dt. 


n=-N 


(9) 


Taking into account that 


/ it) = j fit) + J fi-t) 
0 0 
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the equation 0 can be rewritten 


as 


F(u) 




N 


—27rii't 


dt 


C-v/TT 


n=-N 


+ -^ ^ e-W‘f(-nh) / 


AT 


(10) 


Since change of the variable t/c ^ t leads to 

CX) 

f { ( nht\ 

/ exp I —- j exp I 2^- J exp {±2mth') dt 

0 

OD 

= c y exp (—exp ^2-^ exp (±27ricz/t) dt, 


the equation (10) can be rearranged as 


F(u) 


h 


N 

E- 

n=-N 


“)'/ (nh) / 


e-(-)V(-nh)/e 

n=-w 

Comparing this equation with equation ([^ yields 




F (u) e ^) f (nh) w {—TTciy,—nh/c) 


N 


n=—N 
N 


+ — e ^ ) f {—nh) w {ttcu,— nh/c). 


( 11 ) 


n=-N 


_/ nh ^ 

Lastly, dehning the constants an = he 1 ^ / [nh) /2 we can rewrite approx¬ 

imation 0 in a more compact form as a weighted sum 

^■(0 = -^ {/(()}(>') 

( 12 ) 


n=Af 


[anW {—'Kcv, —nh/c) + a-nW {ttcu, —nh/c)]. 


n=—N 






The equation (12) can be expressed through any functions that have 


been considered in the section 2 above. For example, using the identity ([^ 
after trivial rearrangements we get the Fourier transform in terms of error 
functions of complex argument 

n=N 


F(u) 


n=-N 


~ [1 + erf {nh/c - iTTCu)] 


[1 + erf {nh/c + ivrcz/)] \. 


It can also be shown that the approach based on a weighted sum of the 
complex error functions can be generalized to the Laplace transform. 

The derivation of approximation for the inverse Fourier transform is 
straightforward now. Comparing integrals 0 , (§ and using approximation 
(12) we immediately obtain 


f (t) = {F (u)} (t) 

n=N 

~ (vrct, —nh/c) + (—TTct, 


-nh/c)], 


(13) 


n=-N 


_/ nh “X ^ 

where the coefficients are calculated as a* = he \ ) F {nh) /2. 

Since the Gaussian function ([^ rapidly decreases with increasing |f|, only 
few terms with negative index n actually contribute to shape the curve / {t) 
along the positive t-axis where the Fourier integration takes place accord¬ 


ing to equation (10). As a result, it is sufficient to take into considera 

l2 


tion only, say, hrst three terms he [(*+3^)/d j (cy^), he [h+2/*)/c] j (cy^) and 
}i^-{h+h)ic\ j with negative index n. Thus, the number of the sampling 

points is reduced by almost factor of 2 from 2N -|- 1 to -|- 4. Consequently, 


the equations (12) and (13) can be simplihed as 


F{v)=F{f{t)}{v) 


n=N 


[anW {—'ncu, —nh/c) -|- (ttcz/, —nh/c)] 


(14) 


n=—3 


and 


f {t) = {F {u)] {t) 

n=N 

ft; [a^w {net, —nh/c) + a'_^w {—TTCt, —nh/c)y 


(15) 


n=—3 
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It should be noted that a real-time computation of the complex error 
functions in equations (14) and (15) may not be fast enough for the wavelets 
with extended effective length. However, this problem can be easily resolved 
by using the following precomputed values 

Pn,! = w{-7icu,-nh/c), 

/3n,2 = w{TTCP,-nh/c) 


and 


n* h 

Pn,l = 2 ® ^ ' 

h 


w {iTct, —nh/c), 


Pn,2 = 2 ® W {-7ict,-nh/c) 

that can be stored as datasets in a computer memory. In particular, the 
application of these ready-made precomputed values enables us to simplify 
considerably the approximations (14) and (15) for the forward and inverse 
Fourier transforms as given by 


N 


N 


F{u)=F {/ (t)} (i/) ^ ^ / (nh) /3n,2 (16) 


n=—3 


n=—3 


and 


N 


N 


f (i) = r-' {F (,.)} ^ F (nh) 131, + K .2 , (17) 

n=—3 n=—3 

accordingly. 

3.2 Damping harmonic series 

As an application of the precomputed values (3n,i, /dn ,2 and (3^^, (3^2 ap¬ 


proximations (16) and (17) helps us to avoid a real-time computation of the 


complex error functions, this technique signihcantly accelerates an algorith¬ 
mic implementation. However, there is an alternative way for rapid perfor¬ 
mance in the numerical forward and inverse Fourier transforms. Specihcally, 
using a remarkable property of the complex error function 


w (z 


= 2e- 


w I 


(18) 
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the approximations (12) and (13) can be significantly simplified and ex¬ 
pressed in terms of a damping harmonic series. 

Let us rewrite equation (12) as 

F (z/) [aQW (— TTCz/, 0) -|- UqU) (ttcz/, 0)] 

N 


-|- [dnW {—Tich', —nh/c) + a-nW ( ttcz /, —nh/c)] 

n=l 

-1 

+ [<^nW {—TTCU, —nh/c) + a-nW {'KCl', —nh/c)] 


n=-N 


or 


F (z/) [a^w (—TTCZ/, 0) -f a^w (ttcz/, 0)] 


N 


F [o-nW {-TTCu, —nh/c) + a-nW {ttcu, —nh/c)] 


n=l 

N 


(19) 


(— TTCZ/, nh/c) + anW (ttcz/, nh/c)]. 


71=1 


The identity (18) can be expressed through arguments x and y as 
w {x, y) = — w (—x, —y) = — w (—x, —y). 

This leads to 

w (iTTcz/, -nh/c) = - w (^vrcz/, nh/c). 

Therefore, we have the following identity 


h 


N 


[anW (—TTCZ/, —nh/c) F a-nW (ttcz/, —nh/c)] 


71=1 


= h 


N 


^—{'KCv)^ — 2'Kiunh y ^_ Tlh^ ^—{7vci')^-\-27viunh 


71=1 


N 


— — [anW ( ttcz /, nh/c) F (y.-nW {—ircu, nh/c )], 


71=1 
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/ nn \ 

where hf [nh) = 2Q;„e'' ! . Substituting this identity into approximation 

(|T^ yields 


F (z/) ~ 2 ( —TTCZ^, 0) + w {ttcp, 0)] 

N 


+ [f + / i-nh) . 


n=l 


Since (see [I3] for details) 

w {—X, 0) + tn {x, 0) = 26“"^^, 
this approximation can be expressed as 
F (z/) ^hf (0) 

N 

+ [f + / i-nh) 

n=l 


or 


N 


F(z/) = ^{/(t)}(z/)^he-("-) Y^f{nh)e 


—2Fii/nh 


( 20 ) 


n=-N 


The inverse Fourier transform can be found from approximation (13) by 
using the same derivation procedure 


N 


f (t) = F~^ {F (z/)} (t) ^ ^ 


27 zitnh 


( 21 ) 


n=-N 


Although the approximations (20) and (21) consist of harmonic func¬ 
tions = cos (27rz/nh) — i sin {2'n'unh) and = cos (2Titnh) -|- 

i sin {27rtnh), the resultant forward and inverse Fourier transformed wavelets 
remain, nevertheless, non-periodic due to presence of the damping exponen¬ 
tial functions and . Consequently, unlike the conventional 

discrete Fourier transform the approximations (20) and (21) have no distor¬ 
tion that may appear as a result of aliasing. 
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4 Sample of computation 

If a function / (t) can be separated into the even and odd parts as 


/(f) = /+ (f) + /- (f), 


odd 


where /■*■ (—t) = (t) and / (—t) = — / (f), then according to definition 

([^ the forward Fourier transforms are given by 




oo 

J /■*■ (t) cos [2711^1) dt 

— OO 


( 22 ) 


and 


-^{/ (*)}(") = -' f (t)sm{ 27 rvt)dt. 


(23) 


respectively. Therefore, using approximations ( |20[ ) and (21) for the forward 
Fourier transforms over the even f~^ (t) and odd f~ (t) parts we obtain 


N 


G {u) ^ he ^ f~^ {nh) cos {2Tivnh) 


and 


n=-N 


N 


H (u) ^—ihe / (nh) sin (27rz/nh), 

respectively. Since 


n=-N 


-1 


N 


(nh) cos (27rnnh) = /’*' (nh) cos (27rnnh), 

n=—N n=l 

f~^ (Oh) COS (27rz/0h) = (0) 


and 


Af 


/ (nh) sin (27rnnh) = / (nh) sin (27rz/nh), 

n=—N n=l 

f~ (Oh) sin (27rz/0h) = 0 


(24) 


(25) 
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from equations (24) and (25) it follows that 


and 


G(u)=T{f* (i)}M 

I /■*■ (0) /2 + ^ {nh) cos {^nunh) 


N 


n=l 


= W}M 

N 

~ f~ {nh) sin (27rz/nh). 

n=l 


(26) 


(27) 


fit) 



Fig. 3. The function / (t) defined according to equation (28). 


Consider the following example 

2f + l, 


fit) = 


0 , 


- 1/2 < t < 1/2 

otherwise. 


(28) 


The equation (|28|) can be separated into the even 

1 , 


r it) = 


0 , 


- 1/2 < t < 1/2 

otherwise 


(29) 
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and odd 


r it) 


2t, -l/2<t<l/2 

0, otherwise 


(30) 


parts. Figure 3 shows that the function / (t) while Fig. 4 depicts its even 
/■•■ (t) and odd f~ (t) parts by blue and red curves, respectively. 

From Fig. 4 we can see that the functions {t) and f~ {t) are the 
wavelets (or the pulses) localized within the region t G [—1/2,1/2]. Beyond 
this region the functions {t) and f~ {t) are equal to zero. Therefore, the 
effective length of these wavelets isAf = 2xl/2 = l. Since according to 
approximation ([^ we applied 2iV+l sampling points, the step h between two 
adjacent sampling points can be determined from the formula {2N + 1) h = 
At. Thus, by choosing iV = 50 and N = 300 we can hnd the corresponding 
steps to be h = 0.0099 and h = 0.00166389, respectively. 


f-(t) 



Fig. 4. The even /+ [t] (blue) and odd / [t] (red) parts of the 
function / {t). 


As the function f~ (t) is odd, its Fourier transform H (u) is purely imag¬ 
inary according to equation (23). Figure 5 depicts G (u) and Im [H (z/)] cor¬ 
responding to the functions /+ (t) and f~ (t) by blue and red curves, re¬ 
spectively. The curves for the functions G (z/) and Im [H (z/)] are obtained 
numerically by using equations (26) and (27) at iV = 50, h = c = 0.0099. 

The Fourier transforms for the even f~^ [t) and odd f~ {t) parts of the 
function / (t) can be found analytically. In particular, substituting the equa- 
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G(v), Im[H{v)] 



Fig. 5. The Fourier transform functions corresponding to the even 
(blue) and odd (red) parts of the equation (28). 


tions (29), (30) into approximations (22), (23) and considering the fact that 
these wavelets are not zero-valued only at t G [—1/2,1/2] we get 


and 



1/2 

r 

G{u) 

= / COS 


J 

-1/2 

1/2 

/* 


/ 

2 t sin { 2 ttv, 

-1/2 



TTU COS TTZ/ — Sin TTI/ 


^ 2^2 


Therefore, it is convenient to dehne the differences by using these functions 

N 


Apte = G (u) — 2he I (g) /2 + /'*' {nh) cos {2Tivnh) 

V n=l / 

= sine {tiv) — ( /+ ( 0 ) /2 -|- ^ (nh) cos {2T^vnh) 


n=l 
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and 


N 


Aim = Im [H (z/)] — < —2he / (nh) sin (^Tivnh) 


n=l 


Tiu COS (ttz/) — sin (ttz/) 


N 




— < —2he / {nh) sin {2'Kvnh) 


n=l 


ARe(v) , A,^(V) 



Fig. 6. The differences A^e (blue curve) and Aim (red curve) 
computed at A = 50, = c = 0.0099, respectively. 


Fignre 6 illustrates the differences A^e (blue curve) and Aim (red curve) 
computed at A = 50, h = c = 0.0099. As we can see from Fig. 6, the 
differences A^e and Aim are within the range ±0.001. Further increase of 
the integer N significantly improves the accuracy. This can be seen from 
Fig. 7 showing that the differences Ars (blue curve) and Aim (red curve) 
computed at A = 300, h = c = 0.00166389 remain within the narrow range 
±0.00003. 

5 Conclusion 

We present a new approach for numerical computation of the Fourier integrals 
based on a sampling with the Gaussian function of kind / {c^/^i). It 
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ARe(v), A,^{V) 



Fig. 7. The differences A^e (blue curve) and Aim (red curve) 
computed at A = 300, h = c = 0.00166389, respectively. 


is shown that the Fourier transform can be expressed as a weighted sum of 
the complex error functions. Applying a remarkable property of the com¬ 
plex error function, the weighted sum of the complex error functions can be 
signihcantly simplified as a damping harmonic series. Unlike the conven¬ 
tional discrete Fourier transform, this methodology results in a non-periodic 
wavelet approximation. Therefore, the proposed approach may be practically 
convenient and advantageous in algorithmic implementation. 
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